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Microfluidic  Injector  Models  Based  on 
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Abstract — Lab-on-a-chip  (LoC)  systems  can  be  functionally  decom¬ 
posed  into  their  basic  operating  devices.  Common  devices  are  mixers, 
reactors,  injectors,  and  separators.  In  this  paper,  the  injector  device  is 
modeled  using  artificial  neural  networks  (NNs)  trained  with  finite  element 
simulations  of  the  underlying  mass  transport  partial  differential  equations 
(PDEs).  This  technique  is  used  to  map  the  injector  behavior  into  a  set  of 
analytical  performance  functions  parameterized  by  the  system’s  physical 
variables.  The  injector  examples  shown  are  the  cross,  the  double-tee,  and 
the  gated-cross.  The  results  are  four  orders  of  magnitude  faster  than 
numerical  simulation  and  accurate  with  mean  square  errors  (MSEs)  on 
the  order  of  10-4.  The  resulting  NN  training  data  compare  favorably  with 
experimental  data  from  a  gated-cross  injector  found  in  the  literature. 

Index  Terms — Electrokinetic,  injector,  lab-on-a-chip  (LoC),  microflu¬ 
idic,  neural  network,  simulation. 

I.  Introduction 

Microfluidic  lab-on-a-chip  (LoC)  systems  have  been  studied  for 
more  than  a  decade  and  have  many  applications  in  biology,  medicine, 
and  chemistry  [1],  [2].  They  generally  perform  chemical  analysis 
involving  sample  preparation,  mixing,  reaction,  injection,  separation 
analysis,  and  detection.  Compared  to  traditional  analytical  labs,  LoC 
has  the  significant  advantage  of  increased  analysis  speed,  parallel 
processing,  and  high  integration  and  automation. 

The  simulation  of  complex  LoC  systems  requires  computationally 
expensive  numerical  solutions  to  partial  differential  equations  (PDEs). 
As  the  design  of  LoC  systems  requires  many  repeated  simulations, 
iterative  design  using  numerical  simulation  is  computationally  infea¬ 
sible.  A  much  more  efficient  alternative  involves  functional  decompo¬ 
sition  into  a  series  of  interconnected  blocks,  as  previously  proposed 
for  the  mixer,  injector,  and  separator  [3]— [5],  which  can  be  used 
to  compose  an  entire  LoC  [6].  Behavioral  artificial  neural  network 
(hereafter  referred  to  as  simply  neural  network,  or  NN)  modeling 
makes  such  an  efficient  decomposition  possible.  For  the  mixer  and 
separator,  band  shape  assumptions  are  used  with  analytical  techniques 
to  simplify  the  PDEs  into  several  ordinary  differential  equations. 
For  the  injector,  the  desired  output  quantity  is  a  scalar,  and  not  a 
vector,  as  in  the  mixer  and  separator;  hence,  we  use  numerical  NN 
techniques  instead. 

A  diagram  representing  a  typical  LoC  consisting  of  a  mixer,  reactor, 
injector,  and  separation  system  is  shown  in  Fig.  1.  There  are  two  stages 
to  its  operation.  In  the  first  stage,  a  voltage  is  applied  to  nodes  1,  2, 
3,  4,  and  5.  The  resulting  electric  fields  induce  electrokinetic  flow 
[7],  which  causes  the  fluid  to  be  pumped  through  the  channels.  This 
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Buffer  Sample 


Fig.  1.  Diagram  of  a  canonical  LoC. 

dilutes  the  sample  from  node  2  with  buffer  from  node  1  in  the  mixer 
before  it  is  reacted  with  a  catalyst  and  ultimately  focused  by  nodes 
4  and  5.  In  the  second  stage,  called  electrophoretic  separation,  new 
voltages  are  applied  to  nodes  1-5,  which  causes  a  band  of  the  reacted 
analyte  to  be  injected  into  the  separation  system.  Because  the  analyte 
can  be  composed  of  multiple  species  of  different  charge,  the  individual 
species  migrate  at  different  speeds  in  the  electric  field  causing  them  to 
separate  from  each  other. 

The  importance  of  the  injector  as  a  component  in  a  microfluidic 
separation  system  derives  from  the  fact  that  it  defines  the  shape  and 
quantity  of  analyte  that  will  be  used  for  separation  and  analysis. 
Various  forms  of  microfluidic  electrokinetic  injectors,  such  as  the  tee, 
the  double-tee,  the  cross,  the  double-cross,  and  the  gated-cross,  have 
been  introduced  [8]— [12].  A  first  generation  injector  model  produced 
by  the  authors  was  specific  to  the  cross  injector  and  was  defined  by 
a  two-dimensional  (2-D)  parameter  space  [4].  This  work  improves 
on  previous  generation  injector  modeling  by  using  NN  behavioral 
modeling  concepts  to  create  a  library  of  models  for  injectors,  including 
the  cross,  the  double-tee,  and  the  gated-cross,  each  defined  by  a  four¬ 
dimensional  (4-D)  parameter  space.  The  methodology  is  not  specific 
to  any  one  injector  topology  and  has  been  used  in  macromodels  that 
have  many  more  dimensions  [13]. 

The  methodology  described  here  involves  an  exploration  of  a  rele¬ 
vant  portion  of  the  injector  physical  parameter  space  using  numerical 
solutions  of  the  convection  diffusion  equation.  These  solutions  are  then 
used  to  train  a  NN  that  describes  the  performance  of  the  injector.  The 
behavioral  output  parameters  are  sufficient  to  construct  a  Gaussian 
distribution,  approximating  the  width-averaged  plug  concentrations 
as  a  function  of  longitudinal  position,  for  the  input  to  the  separation 
channel.  Once  the  NN  has  defined  the  behavioral  model,  an  explicit 
analytic  equation  defining  the  NN  is  created  and  can  easily  be  ported 
into  any  modeling  environment.  This  modeling  method  is  related  to  the 
work  done  by  [14],  however,  in  our  method,  we  have  not  utilized  quasi¬ 
random  training  sequences,  because  our  training  data  sets  are  not  large 
enough  to  achieve  the  asymptotic  low-discrepancy  distribution  [15]. 
Furthermore,  our  system  has  a  larger  number  of  physical  variables 
making  the  computational  cost  of  our  network  training  greater,  thus, 
we  use  the  Buckingham-7r  theorem  to  reduce  the  size  of  the  variable 
space  needed  for  training  simulations  [16]. 

The  scope  of  this  work  is  primarily  focused  on  the  application  of 
composable  system  simulation  for  continuous  flow  microfluidic  LoC 
systems.  The  broad  range  of  macromodeling  application  of  NNs  from 
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Fig.  2.  Example  injections  and  geometries  for  (a)  the  cross,  (b)  the  double-tee,  and  (c)  the  gated-cross  injectors.  The  light  shades  indicate  regions  of  high  analyte 
concentration. 


very  large  scale  integrated  (VLSI)  computer-aided  design  (CAD)  [13], 
microfluidics  [14],  and  radio  frequency  (RF)  microelectromechanical 
system  (MEMS)  [17]  suggests  that  further  application  in  microfluidics, 
such  as  digital  microfluidics,  should  be  possible,  but  the  exploration  of 
such  possibilities  are  beyond  the  scope  of  this  work. 

The  injector  topologies  will  be  introduced  in  the  next  section. 
Section  III  will  introduce  the  method  used  to  model  the  injector 
device.  The  physical  variables  of  the  injectors  and  their  nondimen- 
sional  reduction  will  be  shown  in  Sections  III- A  and  III-B,  and  the 
parameters  governing  the  creation  of  the  NN  model  will  be  shown  in 
Sections  III-C  and  III-D.  Finally,  the  impact  of  the  resulting  model  on 
a  composable  system  simulation  will  be  shown  in  Section  III-E,  with 
conclusions  in  Section  IV. 

II.  Injector  Topologies 

The  three  injectors  considered  in  this  paper  are  as  follows: 

1)  the  cross;  2)  the  double-tee;  and  3)  the  gated-cross.  Each  injector 
is  defined  by  its  geometry  and  operating  procedure.  The  geometries 
and  operating  steps  are  shown  in  Fig.  2. 

The  cross  injector  [Fig.  2(a)]  is  a  simple  intersection  of  two  channels 
of  the  same  width.  The  injector  is  operated  by  first  loading  the  sample 
into  the  injection  chamber,  followed  by  a  change  in  voltage  that 
dispenses  a  portion  of  the  injected  flow  into  the  separation  channel.  In 
both  the  loading  and  dispensing  stages,  accessory  fields  can  be  used 
to  manipulate  the  shape  of  the  injected  plug.  In  the  loading  stage, 
accessory  pinching  (from  E2l  and  EAL)  can  be  used  to  achieve  a 
more  compact  plug  and  to  reduce  the  electrokinetic  bias  [8].  In  the 
dispensing  stage,  accessory  pullback  (from  E1D  and  E3D)  prevents 
leakage  into  the  separation  channel. 

The  double-tee  injector  [Fig.  2(b)]  is  a  more  general  form  of  the 
cross  injector  where  the  loading  channels  are  offset  to  produce  a  larger 
injection  plug.  The  operation  scheme  is  the  same  as  with  the  cross,  with 
two  stages  of  operation  and  optional  accessory  fields  for  plug  shaping. 

The  gated-cross  injectors  [Fig.  2(c)]  has  the  same  geometry  as  the 
cross  but  uses  a  three-stage  operating  scheme.  The  gated-cross  control 
scheme  allows  for  continuous  flow  injection  so  that  new  sample  can 
be  loaded  at  the  same  time  that  previously  dispensed  plugs  are  run 
through  the  separation  channel.  The  first  step  of  operation  involves 


creating  the  gate  in  the  injection  chamber  by  counterflowing  an  analyte 
stream  against  a  buffer  stream.  The  second  step  involves  removing 
the  applied  potential  from  the  buffer  port,  which  allows  a  portion 
of  analyte  to  overflow  the  injection  chamber.  The  final  step  is  a 
return  to  the  applied  potentials  of  the  first  stage,  which  reestablishes 
the  gate,  while  simultaneously  injecting  the  overflown  analyte  into 
the  separation  channel.  Since  the  action  taken  in  the  second  stage 
involves  floating  only  one  node  while  all  others  remain  unchanged, 
no  independent  parameters  are  introduced  by  this  stage.  The  fields  and 
flow  patterns  are  completely  determined  by  the  parameters  set  in  the 
first  stage. 

III.  Modeling  Methodology 

The  complexity  of  the  geometry  and  the  resulting  transitional 
electric  field  structure  [18]  makes  behavioral  models  of  the  detailed 
physical  performance  of  the  injectors  very  difficult.  The  modeling 
approach  used  in  this  paper  avoids  the  complexity  of  the  underlying 
physical  structure  while  maintaining  an  accurate  description  of  the 
injector  performance  based  upon  the  key  characteristics  that  define  the 
input/output  mapping. 

The  injector  modeling  methodology  consists  of  four  steps.  1)  The 
Buckingham-7r  theorem  is  used  to  reduce  the  order  of  the  space  of 
physical  parameters  describing  the  behavior  of  the  injector  to  a  reduced 
set  of  nondimensional  parameters.  In  all  three  of  the  injectors  shown 
in  this  paper,  the  number  of  nondimensional  parameters  describing 
the  injector’s  operation  is  reduced  to  four,  which  is  relatively  small 
compared  to  the  dimensions  of  typical  VLSI  analog  design  spaces  [13]. 

2)  Numerical  simulations  are  carried  out  at  points  in  the  nondi¬ 
mensional  parameter  space.  The  goal  is  to  use  as  few  numerical 
simulations  as  possible,  since  they  are  computationally  expensive. 

3)  An  NN  is  constructed  to  analytically  describe  the  parameter 
space.  Functions  available  in  MATLAB  7  [19]  were  used  to  train 
the  NN.  A  feed-forward  back-propagation  NN  topology  was  chosen 
for  its  strong  nonlinear  regression  capabilities  and  its  ability  to  be 
converted  to  an  explicit  algebraic  function.  The  NN  as  a  regressor  has 
the  ability  to  handle  high-dimensional  spaces  with  relatively  sparse 
data  sets  [20].  The  NN  also  “learns”  the  functional  mapping  without 
any  knowledge  of  the  underlying  physics  or  basis  functions.  4)  The 
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TABLE  I 

Summary  of  Nondimensional  Parameters:  The  Top  Four 
Columns  are  Dynamic  and  the  Bottom  Two  are  Geometric 
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trained  network  is  converted  to  an  explicit  algebraic  function,  suitable 
for  use  in  any  software  environment  for  simulation  or  synthesis.  This 
conversion  is  done  with  a  straightforward  parsing  algorithm,  which 
can  be  applied  to  any  generalized  feed-forward  topology.  Thus,  the 
main  contribution  of  this  paper  is  a  method  for  creating  microfluidic 
device  models  using  a  reduced  set  of  nondimensional  variables  and 
input/output  interface  parameters. 


A.  Defining  the  Ti-Space 

Implementation  of  the  Buckingham- tt  theorem  to  the  cross  model 
is  described  in  [4].  Extension  of  the  results  to  the  three  injectors 
described  here  is  summarized  in  Table  I. 

For  all  three  injector  topologies,  the  dynamic  parameters,  7Ti(c,d,g) 
through  7t4(c,d,g)>  are  the  independent  variables  used  to  create 
the  models.  The  geometric  parameters  are  held  constant,  such  that 
tt5(c, d,g)  =  1  and  7t6D  =  2.  These  geometries  are  representative  of 
the  vast  majority  of  designs  fabricated  to  date  for  each  type. 

Physically,  the  cross  and  the  double-tee  are  very  similar.  From 
a  synthesis  point  of  view,  the  cross  is  simply  a  special  case  of 
the  double-tee,  where  7t6d  =  0.  For  both  of  these  topologies,  the 
7Ti(c,d)  and  7t2(c,d)  parameters  describe  the  ratio  of  the  accessory 
fields  to  the  driving  fields.  For  the  loading  stage,  this  describes  the 
amount  of  pinching  that  is  applied  to  the  incoming  analyte  stream,  and 
for  the  dispensing  stage,  this  describes  the  amount  of  pullback  applied 
to  the  dispensed  band  [4].  The  7t3(c,d)  and  7t4(c,d)  parameters  de¬ 
scribe  the  Peclet  number  for  each  stage. 

The  gated-cross  has  a  set  of  parameters  that  differ  from  those  of 
the  cross  and  the  double-tee.  The  first  parameter  n1G  represents  the 
extent  to  which  the  gate  is  closed.  As  discussed  in  [12],  as  long  as 
E\  >  E2,  the  gate  will  remain  closed  in  the  limit  of  no  diffusion.  As 
7Tig  is  reduced,  the  gate  is  further  closed,  so  in  practice,  Ei  >  E2. 
The  second  parameter  represents  the  ratio  of  the  buffer  electric  fields 
(E3,EA)  to  the  analyte  electric  fields  (E2,Ei).  As  7t2g  increases, 
the  buffer  fields  become  larger  relative  to  the  analyte  fields.  The  third 
parameter  7t3g  represents  the  Peclet  number  of  the  loading  phase.  The 
final  dynamic  parameter  7t4g  is  a  measure  of  the  ratio  of  length  of 
the  injected  plug,  as  related  to  the  floating  time  TLD,  to  the  channel 
width. 

A  significant  concern  in  the  operation  of  a  gated-cross  injector  is 
the  leakage  that  can  occur  in  the  separation  channel  if  the  gate  is  not 
sufficiently  closed,  as  seen  in  Fig.  3.  If  the  leakage  is  too  great,  the 
injector  will  not  operate,  because  the  increased  noise  floor  will  make 
a  separation  impossible.  A  region  of  feasibility  must  be  determined 
within  which  the  injector  can  be  modeled.  This  leakage  was  analyzed 
in  [12]  as  a  function  of  the  gate  closure  and  the  system  Peclet  number. 
They  determined  a  boundary  indicated  by  1%  leakage  of  the  flux  of  the 
analyte  from  the  source  reservoir  into  the  separation  channel.  In  this 
work,  we  determine  the  boundary  for  a  more  complete  set  of  physical 
parameters  to  create  a  region  of  feasibility. 


Fig.  3.  Gated-cross  leakage  tested  at  Peclet  numbers  from  19  at  contour 
1  to  oo  at  contour  5.  The  contours  are  defined  by  the  7%  of  the  maximum 
concentration  and  agree  very  well  with  the  numerical  and  experimental  results 
found  in  [12]. 

B.  Simulating  the  Ti-Space 

The  injectors  were  simulated  in  FEMLAB  using  the  convection 
diffusion  equation  [21]  given  as 

dc  —> 

—  +u-Vc  =  ftV2c  (1) 

at 

where  c  is  the  concentration  distribution,  u  is  the  velocity  field,  and  k 
is  the  diffusion  coefficient. 

To  solve  for  the  electrokinetic  flow,  first,  Laplace’s  equation  is 
solved  for  the  electric  potential  as 


V2<f>  =  0  (2) 

which  is  used  to  determine  the  velocity,  which  is  proportional  to  the 
electric  field,  through 


u  =  — /iVT>  (3) 

where  p  is  the  electrokinetic  mobility,  and  the  negative  gradient  of  the 
electric  potential,  T>,  is  the  electric  field. 

In  all  simulations,  the  channels  are  long  enough  to  contain  the 
transitions  in  the  electric  field  structure.  The  buffer,  analyte,  and  waste 
channels  were  all  of  the  same  length,  4if,  for  each  injector  topology. 
The  separation  channel  length,  set  so  that  it  can  encompass  the  entire 
injected  band,  is  9 w  for  the  cross,  12w  for  the  double-tee,  and  13if 
for  the  gated-cross.  This  will  allow  the  numerical  simulation  data  to 
be  used  to  train  the  NN  to  contain  a  highly  accurate  description  of  the 
transitional  electric  field  structure  and  its  effect  on  the  dispersion  of 
the  analyte. 

For  each  injector  topology,  using  conservation  of  current  and  as¬ 
suming  that  all  channel  widths  are  the  same,  a  relationship  between 
the  steady  electric  fields  can  be  determined.  As  in  prior  work  on  the 
cross  injector  [4],  the  electric  field  magnitudes  in  the  cross  and  the 
double-tee  are  related  for  the  loading  stage  by 

E3L  =  E1L  +  2  E2,4l  (4) 

and  for  the  dispensing  stage  by 

E^d  —  E2d  +  2  Ei^d  (5) 

where  it  is  assumed  that  the  accessory  fields  are  applied  symmetrically, 
i.e.,  £2L  =  £4L  =  E2,4l  and  E±r>  =  E3r>  =  Elj3r>.  The  fields  for 
the  gated-cross  are  related  by 


E2  +  E3  —  EA  -f-  E\ . 


(6) 
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TABLE  II 

7t-Space  Domains  for  the  Training  Simulations, 
Covering  a  Wide  Range  of  Physically  Viable  Values 


7lX 

^2 

-73 

|. 

Cross 

[0:1/2] 

[0:2] 

[50:5000] 

[50:5000] 

Double-Tee 

[0:1/41 

[0:41 

[50:50001 

[50:50001 

Gated-Cross 

[1/10:11 

[1:4] 

[50:50001 

[1/20:1/31 

Using  (4)-(6),  and  the  conversion  matrix  introduced  in  [4],  it  is 
possible  to  calculate  the  electric  potentials  for  the  boundary  conditions, 
as  they  change  throughout  the  7r-space. 

To  perform,  on  average,  200-300  simulations  per  injector,  a 
FEMLAB  script  was  written  to  automate  the  setup  and  execution  of 
the  simulations.  The  7r- space  domains  for  each  of  the  injectors  are 
defined  in  Table  II.  A  parallel  computing  cluster  was  utilized,  which 
reduced  several  weeks  of  computation  to  several  days  on  a  shared 
100-node  Beowulf  cluster.  This  is  a  one-time  computational  expense 
used  to  create  the  much  more  efficient  NN  description,  which  “learns” 
from  the  results  of  the  simulations. 

The  FEMLAB  algorithms  have  been  validated  against  many  ex¬ 
periments  found  in  the  literature  demonstrating  microfluidic  mixing, 
joule  heating,  injection,  and  separation  [3],  [5],  [22].  To  validate  the 
simulations  for  the  new  gated-cross  topology,  a  test  was  run  to  measure 
the  leakage  of  the  injector  during  the  loading  stage,  as  was  first  done  in 
[12].  Fig.  3  shows  the  results  of  the  simulations  at  the  same  parameters 
found  in  [12].  The  resulting  isoconcentration  contours  agree  very  well 
with  [12]  numerical  and  experimental  results. 

After  numerically  simulating  points  throughout  the  parameter  space, 
the  form  of  the  output  of  the  behavioral  model  is  chosen.  Analysis 
of  the  numerical  results  can  be  done  to  extract  any  desired  behavioral 
outputs.  In  this  work,  the  outputs  were  chosen  to  be  the  peak  height  and 
variance  of  an  effective  Gaussian  describing  the  transversely  averaged 
concentration  profile  output  from  the  injector,  such  that  the  area  and 
variance  of  the  original  distribution  are  conserved.  The  resulting  form 
with  nondimensional  outputs  is 


/  (^l(fc),^2(fc)5^3(/e)?^4(/e)) 


geff 


L  C0 


(k) 


(7) 


where  k  E  [C,D,G\,  and  CQ  is  the  value  of  the  uniform  concentration 
of  the  input  species  from  the  reactor. 

While  the  NN  is  inherently  capable  of  describing  multiple-input 
multiple-output  models,  we  use  a  multiple-input  single-output  NN  to 
train  each  output  separately. 

Using  effective  Gaussian  outputs  as  the  interface  parameters  for  a 
system  simulation  accomplishes  two  purposes,  namely:  1)  the  system 
simulation  is  efficient;  and  2)  while  there  is  sufficient  information  to 
accurately  describe  the  band  after  it  has  traveled  a  short  distance  down 
the  separation  channel,  there  is  not  enough  information  to  describe  the 
non-Gaussian  shape  of  the  band  with  high  accuracy  immediately  after 
it  leaves  the  injector. 

As  long  as  diffusion  transversely  homogenizes  the  band  before 
detection  occurs,  or  before  another  source  of  dispersion,  such  as  a 
channel  bend,  is  introduced,  the  approximation  will  provide  accurate 
results.  By  requiring  the  time  it  takes  particles  to  diffuse  across  the 
channel  to  be  shorter  than  the  time  it  takes  the  particles  to  reach  the 
first  channel  bend  or  detector,  a  bound  is  defined  as 


T  — 


Peu 

L 


<  2 


w 


(8) 


where  Pew  is  the  Peclet  number  based  on  the  channel  width  w  and 
L  is  the  length  of  straight  channel  before  the  detector  or  first  channel 


L  =  0  :  r=  oo  L  =  4.7  mm  :  r=  2 


Fig.  4.  Actual  output  of  the  double-tee  injector  (top  channel)  compared  to 
effective  parameterized  output  of  the  analytical  model  (bottom  channel).  The 
band  on  the  right  traveled  4.7  mm,  which  is  when  r  =  2.  For  these  simulations, 
7Tid  =  186, 7T2D  =  186, 7T3D  =  1/8,  and  7T4d  =  0.57. 


''[output 
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} 


input 
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Fig.  5.  (a)  Activation  of  a  node  is  a  function  of  the  weighted  inputs  of  the 

previous  layer’s  activations,  (b)  General  feed-forward  NN  topology  utilizing 
multiple  hidden  layers  with  an  unspecified  number  of  nodes  with  nonlinear 
activation  functions  and  a  single  output  node  using  a  linear  activation  function. 


bend.  This  expression  is  conservative  in  that  it  requires  the  particles 
travel  across  the  entire  channel  width  to  homogenize  the  band.  In  many 
cases,  as  seen  in  Fig.  2,  the  band  coming  out  of  the  injector  is  axially 
compact  and  fairly  uniform  across  the  channel  and  so  homogenization 
occurs  much  more  rapidly. 

An  example  for  a  double-tee  injection  is  shown  in  Fig.  4.  The  picture 
on  the  left  shows  the  band  immediately  after  the  injector,  and  the 
picture  on  the  right  shows  the  band  after  traveling  4.7  mm  where,  from 
(8),  r  =  2.  As  an  example  of  an  actual  experimental  system,  the  LoC 
developed  by  Chiem  and  Harrison  [23]  has  r  =  0.2,  which  is  an  order 
of  magnitude  within  the  bounds  defined  by  (8). 


C.  Training  the  NN 

An  NN  is  a  mathematical  structure  that  is  adept  at  learning  nonlin¬ 
ear  functional  relations  and  complex  item  categorizations  [20].  The 
general  principle  of  operation  of  the  feed-forward  NN  is  shown  in 
Fig.  5(a),  and  the  topology  of  the  two-layer  NN  used  to  model  the 
injectors  for  this  work  is  shown  in  Fig.  5(b).  The  feed-forward  NN 
is  a  set  of  nodes  that  are  connected  only  to  the  layers  above  and 
below.  Each  node  takes,  as  the  argument  of  an  activation  function, 
the  weighted  sum  of  the  outputs  of  the  layer  below.  To  create  an  NN 
model,  the  number  of  layers,  nodes  per  layer,  node  transfer  function, 
and  interconnecting  weights  must  be  selected. 

It  has  been  shown  that  to  perform  nonlinear  regression  to  an 
arbitrary  degree  of  accuracy,  where  the  number  of  nodes  is  not  an 
issue,  two  layers  are  sufficient  (one  hidden,  one  output)  [24].  Likewise, 
with  classification  networks,  it  has  been  shown  that  a  decision  bound¬ 
ary  of  arbitrary  complexity  can  be  defined  using  only  three  layers 
(two  hidden,  one  output)  [20].  Therefore,  our  regression  models  use 
two-layer  topologies,  and  our  classification  networks  use  three-layer 
topologies. 
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TABLE  III 

Neural  Network  MSE  and  Cross  Validation:  For  All  Three  Injectors,  Both  Validation  Parameters  are  Very  Small 


MSE 

KFCV 

a2/w2 

C 

v-/max/v-'0 

k 

g2/w2 

C  /C 

v-/max/  '-/o 

Cross 

1.2e-4 

2e-5 

10 

2e-6 

9e-5 

Double-Tee 

1.8c-5 

le-5 

8 

3e-6 

2e-5 

Gated-Cross 

4e-5 

4.2e-5 

10 

l.le-4 

2.4e-4 

Fig.  6.  (a)  Variance  and  (b)  effective  peak  concentration  results  for  the  cross,  where  7Tic  and  7T2c  are  fixed  at  0.2.  (c)  Variance  and  (d)  effective  peak  concentration 

results  for  the  double-tee,  where  7Tid  and  7T2d  are  fixed  at  0.05  and  0.3,  respectively.  The  top  plots  show  the  FEMLAB  simulations  at  a  coarse  set  of  points  and 
the  bottom  plots  show  the  much  faster  and  higher  resolution  evaluations  produced  by  the  NN. 


If  the  network  is  to  perform  nonlinear  regression,  the  hidden  layer 
consists  of  bound  nonlinear  activation  functions,  and  the  output  nodes 
are  unbound  linear  activation  functions.  If  the  network  makes  discrete 
classifications,  then  all  nodes  are  bound  nonlinear  functions.  In  both 
cases,  the  actual  shape  of  the  nonlinear  function  is  not  important  as 
long  as  it  is  bound  [20].  The  function  used  in  this  work  is  written  as 


/(*)  = 


2 

1  +  exp(— 2x) 


-  1. 


(9) 


The  number  of  input  nodes  is  determined  by  the  number  of  indepen¬ 
dent  parameters  in  the  model  and  the  number  of  nodes  in  the  output 
layer  is  determined  by  the  number  of  dependent  parameters  in  the 
model.  The  number  of  nodes  in  the  hidden  layer  is  arbitrary.  Nodes 
are  added  or  subtracted  until  the  performance  reaches  a  threshold.  The 
hidden  layers  of  the  injector  NNs  contain  25  nodes. 

After  the  topology  and  initial  node  count  have  been  determined, 
the  network  is  trained  on  the  FEMLAB  simulation  data.  The  process 
of  training  an  NN  involves  modifying  the  weights  between  layers 
until  the  outputs  of  the  network  sufficiently  match  the  simulation 
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Fig.  7.  (a)  Gated-cross  feasibility  space  (dark  is  feasible),  (b)  Variance  and  (c)  effective  peak  concentration  results  for  the  gated-cross,  where  tt^g  and  7T4q  are 

fixed  at  232  and  0.1208,  respectively. 


training  data.  This  takes  the  form  of  an  optimization  problem.  The 
Levenberg-Marquardt  algorithm  included  in  MATLAB  was  used  due 
to  its  speed  for  networks  of  complexity  similar  to  the  ones  shown  here. 
Several  metrics  are  used  to  measure  the  quality  of  the  trained  network. 
To  ensure  that  the  network  is  not  underfit,  the  mean  square  error  (MSE) 
is  minimized  during  training,  and  to  ensure  that  the  network  is  not 
overfit,  i.e.,  has  poor  generalization,  the  trained  network  is  tested  using 
data  that  have  not  been  seen  before.  To  test  for  overfitting  without 
creating  a  separate  validation  data  set,  k- fold  cross  validation  (KFCV) 
is  used  [20]. 

For  the  models  in  this  paper,  no  more  than  350  iterations  were  used 
to  achieve  an  MSE  accuracy  of  less  than  1  x  10-4.  To  improve  the 
network’s  convergence  rate,  the  range  of  the  input  data  is  prescaled. 
In  particular,  the  Peclet  number  is  varied  on  a  logarithmic  rather  than 
linear  scale,  and  the  variance  training  data  for  the  gated-cross  injector 
are  normalized  to  unity.  To  validate  the  NN  models,  first,  the  MSE 
and  KFCV  are  examined;  the  results  are  summarized  in  Table  III,  with 
very  small  validation  parameters  throughout.  As  a  visual  confirmation, 
response  surfaces  are  shown  for  portions  of  the  7r-space  for  each 
injector.  Fig.  6  shows  the  results  for  the  cross  and  the  double-tee 
injector,  and  Fig.  7  shows  the  results  for  the  gated-cross  injector. 

The  graphs  at  the  top  represent  the  sparse  set  of  points  simulated 
in  FEMLAB,  while  those  at  the  bottom  show  the  densely  simulated 
points  using  the  much  faster  NN  after  being  trained  on  the  FEMLAB 
data.  The  NN  operates  in  a  fraction  of  a  second  per  point,  more  than 
four  orders  of  magnitude  faster  than  numerical  simulation. 

Fig.  6(a)  and  (b)  shows  the  variance  and  effective  peak  concentration 
results,  respectively,  for  the  cross,  and  Fig.  6(c)  and  (d)  shows  the 
same  results  for  the  double-tee.  When  the  cross  and  the  double-tee 


are  operated  with  the  same  accessory  fields,  the  cross  has  a  smaller 
variance  and  a  smaller  peak  concentration. 

Fig.  7(a)  shows  the  feasibility  regions  for  the  gated-cross,  and 
Figs.  7(b)  and  6(c)  show  the  variance  and  peak  concentration  results 
for  the  gated-cross  injector,  respectively.  The  feasibility  space  is 
independent  of  7T4g,  because  the  leakage  is  determined  in  the  loading 
stage  independently  of  how  the  band  is  dispensed.  The  most  notable 
feature  of  the  feasibility  space  is  that  as  7t3g,  the  Peclet  number, 
becomes  smaller,  the  amount  of  feasible  space  decreases  significantly. 
This  is  due  to  the  fact  that  regardless  of  how  much  the  gate  is  closed, 
when  diffusion  is  large,  the  diffusive  flux  leaking  into  the  separation 
channel  is  large.  Fig.  7(b)  shows  the  variance  results  in  the  7r-space  of 
the  gated-cross,  and  Fig.  7(c)  shows  the  effective  peak  concentration 
results.  The  missing  surface  areas  in  Fig.  7(b)  and  (c)  represent 
infeasible  regions  as  described  in  Section  III-A. 


D.  Extracting  an  Analytic  Equation 

Since  the  NN  is  connected  in  a  feed-forward  configuration,  it  is 
straightforward  to  identify  the  explicit  relation  for  the  output  node  in 
terms  of  all  the  nodes  in  the  previous  layers.  The  general  explicit  form 
for  the  two-layer  topology  with  a  single  output  node  is 


'  M 

N  1 

Vk  —  Qoxit 

^  ^  Vkj  Qhid 

.3=0 

_i=Q  J_ 

where  y is  the  output  of  the  network,  gout  [£]  is  the  linear  output  node 
activation  function,  ghid[C]  is  the  nonlinear  hidden  node  activation 
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function  from  (9),  and  Vkj  and  Wji  are  the  weights  between  the  hidden- 
output  and  input-hidden  layers,  respectively. 

Once  the  network  is  trained  and  the  weights  are  fixed,  a  sim¬ 
ple  algorithm  can  be  written  to  extract  the  network  data  from  the 
MATLAB  data  structures.  The  result  is  then  used  in  VerilogA  to 
perform  analytical  system  simulations. 


E.  Injector  Model  Impact 

As  a  simple  example  to  illustrate  the  benefit  of  using  the  accu¬ 
rate  NN  model,  a  standard  cross  separation  is  performed  in  Fig.  8. 
Comparison  is  made  between  numerical  simulation  of  the  separation 
using  FEMLAB,  analytical  simulation  using  our  VerilogA  simulator 
with  the  device  models  of  the  injector  and  separation  channel,  and  a 
VerilogA  simulation  using  a  simplified  injector  model.  The  simplified 
injector  model  assumes  that  the  resulting  injection  will  be  a  Gaussian 
approximation  to  a  rectangular  plug  with  concentration  CQ  and  stan¬ 
dard  deviation  equal  to  half  the  channel’s  width,  w/2.  The  Gaussian 
will  have  the  same  standard  deviation  as  the  rectangular  plug  and  a 
peak  height  that  conserves  area,  Cp  —  CQ/  y/ 7t/2 . 

To  compare  the  simulation  results,  the  resolution  of  the  separation 
is  calculated.  Resolution  is  defined  as 


(d2  —  di) 
4crave 


(11) 


where  d2  and  d\  are  the  centers  of  the  two  bands  in  time,  respectively, 
and  crave  is  the  average  standard  deviation  of  the  two  bands  in  time. 
These  quantities  are  meaningful  when  the  electropherogram  peaks  are 
approximated  to  be  Gaussian,  as  they  are  in  this  case.  Four  times  the 
standard  deviation,  therefore,  accounts  for  more  than  95%  of  the  total 
mass  of  the  band.  A  desirable  resolution  is  usually  at  least  1.5,  so  the 
ease  of  the  separation  in  the  example  of  Fig.  8  is  illustrated  by  its  very 
high  resolution  number. 

The  results  show  the  simulation  with  the  NN  injector  model  in 
close  agreement  with  the  numerical  simulation.  The  simulation  using 
the  approximate  injector  model  shows  more  than  25%  difference  in 
resolution.  In  this  system,  the  two  bands  are  still  easily  resolved, 
however,  in  realistic  systems  with  tens  to  hundreds  of  species,  this 
difference  becomes  much  more  significant. 


IV.  Conclusion 

This  work  has  shown  a  methodology  that  was  used  to  model  the 
injector  component  of  a  microfluidic  LoC,  and  examples  were  shown 
for  the  cross,  the  double-tee,  and  the  gated-cross  injector  designs.  The 
key  contributions  of  this  work  were  the  use  of  nondimensionalized 
reduced  variable  sets  for  the  numerical  simulation,  and  the  use  of 
simplified  Gaussian  interface  parameters  to  describe  the  injection  out¬ 
put  to  the  downstream  separation  channel  models.  The  impact  of  the 
NN-based  injector  model  was  seen  to  be  significant  compared  to  a  sim¬ 
ple  approximate  injector  model  and  could  become  a  significant  design 
factor  in  realistic  systems  separating  tens  to  hundreds  of  species.  These 
injector  models  combined  with  other  component  models  allow  for  a 
functional  block  decomposition  of  the  system  for  efficient  simulation. 
The  speed  and  accuracy  of  these  analytic  block  models  present  a  far 
more  feasible  method  of  computer-aided  design  (CAD)  than  using 
numerical  solutions  of  partial  differential  equations  (PDEs)  for  the 
whole  complex  system. 

This  work  lays  the  foundation  for  future  enhancements  that  include 
the  geometric  7r-space  variables  for  a  more  complete  synthesis  formu¬ 
lation.  Due  to  its  ability  to  model  arbitrary  functional  forms  with  an 
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Fig.  8.  Simulations  of  the  cross  separation  system  using  component  sepa¬ 
ration  system  model  and  1)  an  approximate  injector  as  described  in  the  text, 
2)  the  NN  injector,  3)  full  finite-element  method  (FEM)  simulation  of  injector 
and  separation  channel. 


arbitrary  number  of  inputs  and  outputs,  with  further  investigation,  the 
methodology  described  here  can  be  used  to  model  other  microfluidic 
components  for  use  in  composable  system  simulation. 
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